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Abstract 

Molecular dynamics simulations were used to determine the binding affinities between 
the hormone 17/3-estradiol (E2) and different estrogen receptor (ER) isoforms in the rainbow 
trout, Oncorhynchus mykiss. Previous phylogenetic analysis indicates that a whole genome 
duplication prior to the divergence of ray-finned fish led to two distinct ER/3 isoforms, ER/31 
and ER/32, and the recent whole genome duplication in the ancestral salmonid created two 
ERa isoforms, ERal and ERa2. The objective of our computational studies is to provide 
insight into the underlying evolutionary pressures on these isoforms. For the ERa subtype 
our results show that E2 binds preferentially to ERal over ERa2. Tests of lineage specific 
dN/dS ratios indicate that the ligand binding domain of the ERa2 gene is evolving under 
relaxed selection relative to all other ERa genes. Comparison with the highly conserved 
DNA binding domain suggests that ERa2 may be undergoing neofunctionalization possibly 
by binding to another ligand. By contrast, both ER/31 and ER/32 bind similarly to E2 and 
the best fitting model of selection indicates that the ligand binding domain of all ER/3 genes 
are evolving under the same level of purifying selection, comparable to ERal. 



1 Introduction 

Estrogens are essential endogenous hormones that modulate the development and homeostasis 
of a wide range of target tissues, such as the reproductive tracts, breast and skeletal system 
[1]. Estrogenic hormones have multi-faceted and wide-ranging effects in vertebrate animals. 
For estrogens such as 17/3-estradiol (E2) to exert their biological effects, they must interact 
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with cellular estrogen receptors (ER). Studies have shown that ERs are part of two distinct 
estrogenic transduction pathways. One pathway provides a rapid, nongenomic pathway initiated 
by membrane bound ERs at the cell surface IH [21 13] . The other pathway provides direct genomic 
control in which ERs act as transcription factors within the cell nucleus [H [5j . These ERs are 
members of the nuclear receptor superfamily of ligand- modulated transcription factors [H [71 [8] . 
There are two different subtypes of these ERs, referred to as a and /3, each encoded by a 
separate gene. 

Recently, Nagler et al [7J reported the novel ERa2 and both ER/3 isoforms in the rainbow 
trout, Oncorhynchus mykiss, and performed a comprehensive phylogenetic analysis with all 
other known fish ER gene sequences. Their phylogenetic analysis indicates that the duplication 
leading to the two ER/3 isoforms arose prior to the divergence of the ray finned fish attributable 
to a whole genome duplication that occurred in the Teleost ancestor (see Fig |4| ^ . The ERa 
isoforms, on the other hand, appear to have arisen as a result of a second more recent whole 
genome duplication event that occurred in the salmonid ancestor 25-100 million years ago 
|10j . These results indicate that the second ERa isozyme that arose during the earlier genome 
duplication appears to have been lost subsequently, since no other ray finned fish are known to 
have a second ERa isoform. This also indicates that the expected duplications of ER/31 and 
ER/32 were lost subsequent to the salmonid genome duplication. 

The purpose of the study is to employ molecular dynamics simulations to determine the 
binding affinities between E2 and ERs of the different isoforms in the rainbow trout and to use 
the results to provide insight into the underlying evolutionary selection pressure on the ERs. 
Our binding affinity results obtained from insertion and deletion are very similar indicating 
that our simulations are well converged and that accurate estimates of binding affinities were 
obtained. Our results show that E2 binds preferentially to ERal over ERa2. By contrast, 
the difference in binding affinity is less significant for the /3 subtype, i.e., both isoforms bind 
similarly to E2. We also computed dN/dS ratios for the ER isoforms. These results suggest 
that the ERal gene is evolving under relaxed selection compared to all other salmonid ERa 
genes. 

Materials and Methods 

Receptor Structures 

The initial coordinates for the estradiol were first extracted from the human ER-E2 complexes 
(PDB: IQKU) (Fig[l]). Thetopolo gies were then generated by the PRODRG server |llj with the 
options of full charges and no energy minimization. The rainbow trout ER holo structures for 
the E2 binding domain were generated by SWISS-MODEL [12] using human ER as templates 
(PDB entries 1A52 for ERa's and 3ERT for ER/3's). Sequence identities between trout and 
human estrogen binding domains are within the range of 75-85%. The estradiol was first docked 
into the binding pocket of the receptor holo structure with AutoDock [13]. In this protocol, 
the receptor structure is held rigid and the estradiol is free to rotate and explore most probable 
binding poses using the Lamarckian genetic algorithm. The number of genetic algorithm runs 
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Figure 1: Crystal structure of human estrogen receptor binding domain bound to the hormone 
17/3-estradiol. Similar human ER structures were used as templates to generate structures for 
trout ERs. Image was rendered using VMD [1^. 

was set to 1,000 with a population size of 5,000 individuals and 5,000,000 generations. The 
number of evaluations was set to 2,500,000 for each individual in the population to ensure 
thorough exploration of the search space. The mutation rate was set to 0.02 and crossover 
0.8. Two-point crossover was used to generate the offspring at each successive generation. The 
genetic algorithm automatically preserved the 10 best-fit individuals to the next generation and 
the 10 least-fit individuals were not used to generate offspring. A total of 1,000 independent 
docking trials were performed for each of the four ERs. The best binding pose from each trial 
was collected and ranked based on the scores. These best-fit binding poses were first visually 
inspected for consistency with human ER and the one with the highest score was then used as 
the starting structure for the simulations. 

Thermodynamic Cycle 

To estimate E2-ER binding affinities we note that, since the free energy is a state function, it 
permits the selection of an arbitrary path connecting the bound and unbound states. There- 
fore, we decomposed the binding free energy calculation into several steps in which the E2 is 
annihilated (i.e., decoupled) from its bound state in the receptor complex and then made to 
reappear in solution to complete the thermodynamic cycle. For brevity, we subsequently define 
deletion to be when interactions between E2 and its environment are turned off and insertion 
to be when these interactions are turned on. 

Fig [3] shows the thermodynamic cycle we used to calculate binding affinities (see also Refs 
[T5I [TBI El IHl CH]). Starting with upper right schematic and moving clockwise, the fully 
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Figure 2: Restraints used for 17/3-estradiol. We used six harmonic restraints acting on three 
anchor atoms in the E2 and three in the receptor. The restraints consisted of one distance 
(r), two angle (^a and ^b), and three dihedral {4>a, <^B) and ^c) restraints that determine the 
orientation of three carbons in the E2 relative to three a-carbons in the receptors. 

interacting E2 (blue) is first restrained in the binding pocket of the receptor. Here RE represents 
the solvated complex of the E2 and receptor, and AG^^ denotes the free energy of restraining 
the E2 in the binding pocket of the receptor which will depend on the details of restraint. Next, 
the electrostatic and Lennard-Jones interactions of the E2 are gradually turned off (white) 
in two separate steps using alchemical simulations. The free energies AG^^ and AG^f" are 
associated with deleting or inserting electrostatic and Lennard-Jones interactions respectively. 
With the E2 fully decoupled from its environment, the restraint is then removed. The free 
energy AG^g^ is associated with the removal of this restraint. Next the E2 interactions are 
turned back on with no receptor present. The free energy AG^g^. and AGlj are associated with 
turning on the electrostatic and Lennard-Jones interactions respectively. Finally, we account 
for the difference between the standard (Vq) and simulation volume (Kiim)- The binding affinity 
between the estrogen receptors and E2 is thus the sum of the free energies, 

- AGbind = AGfJ + AGi^, + AGff + AGg,, + AG^, + AGf j + AGf„ . (1) 
Restraints 

To facilitate convergence restraints were applied to restrict the positions of E2 relative to the 
receptors. Boresch et al [20j and Mobley et al ^\ reported that the presence of multiple 
metastable ligand orientations can cause convergence problems for free energy estimates. The 
authors further suggested using a restraining potential to keep the ligand in the binding site 
during the simulation process. With such a restraining potential the ligand is no longer required 
to sample the entire simulation volume (particularly a problem when ligand is decoupled). 
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Figure 3: Thermodynamic cycle used for calculating binding affinities. Since the free energy is 
a state function, the calculation of binding affinity is decomposed into several steps [21]. Eqn 
[T]was used to calculate the binding affinity AGbind between the ERs and E2. The gray curved 
rectangle represents the receptor, the blue circle represents E2 with all interactions turned on, 
and the white circle indicates that all interactions are turned off. The spring represents the 
restraints between E2 and receptors. 

Moreover, the restraint minimizes the detrimental effects of end-point singularities commonly 
reported in alchemical simulations |16| [T71 I19j . Mobley et al [T7] also pointed out that the 
equilibrium geometry of the restraints is arbitrary and will not affect the asymptotic estimate 
of the binding free energy. In this work, we judiciously selected anchor atoms from the more rigid 
alpha helices that form the E2 binding pocket. The restraints included one distance (with the 
force constant of 1000 kJ/mol/nm^), two angle (1000 kJ/mol/rad), and three dihedral restraints 
(1000 kJ/mol/rad^) that determine the orientation of three carbons in the E2 relative to three 
a-carbons in the receptors (see Fig [2]). 

Simulation Protocols 

All simulations were performed with the GROMACS 4.0 [22] compiled in single-precision mode 
at a constant temperature of 277 K in a periodic box with an edge length of approximately 8.2 
nm and the default GROMOS-96 43A1 forcefield [23] . The simulation systems each contained 
approximately 16,500 Simple Point Charge (SPC) water molecules [23]. Short-range interactions 
were evaluated using a neighbor list of 1.0 nm updated at every 10 steps. Van der Waals 
interactions used a cutoff with a smoothing function such that the interactions slowly decayed 
to zero between 0.75 nm and 0.90 nm. A long-range analytical dispersion correction was applied 
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to the energy and pressure to account for the truncation of the Lennard- Jones interactions [25] . 
Electrostatic interactions were evaluated using the particle mesh Ewald (PME) |l26i| with a real 
space cutoff of 1.0 nm, a spline order of 6, a Fourier spacing of 0.1 m, and relative tolerance 
between long and short range energies of 10^^. All bonds to hydrogen were constrained with 
LINCS [27] with an order of 12, and a time step of 2 fs was used for dynamics. 

For equilibration, the systems were first minimized using 1,000 steps of L-BFGS (Broyden- 
Fletcher-Goldfarb-Shanno) [28], followed by 1,000 steps of steepest descent minimization. The 
system was then subject to 1.0 ns of simulation using isothermal molecular dynamics. This was 
followed by another 1.0 ns of simulation using isothermal-isobaric dynamics with the Berendsen 
barostat with a time constant of 1.0 ns. For all simulations the temperature was maintained 
at 277 K using Langevin dynamics [29] with a friction coefficient of 1.0 amu/ps. The coupling 
time was set to 0.5 ps, and the isothermal compressibility was set to 4.5 x 10~^ bar~^. 

After equilibration, production simulations were run with isothermal-isobaric conditions 
using Langevin dynamics at the temperature of 277 K. The pressure was maintained at 1.0 atm 
using the Parrinello- Rahman algorithm [30] . The temperature was chosen as it closely resembles 
the water temperature for the natural habitat of rainbow trout. Energies were recorded every 
0.2 ps during production runs, and trajectory snapshots every 1.0 ps. The first 50% of each 
simulation was discarded for equilibration. 

Free Energy Calculations 

We used the formula suggested by Boresch et al |2Dj to analytically calculate the free energy 
AG^st associated with adding the restraints to E2 when decoupled from its environment. We 
also analytically calculated the free energy AGy^ that accounts for the difference between the 
standard (Vq) and simulation volume (Kiim) [20] . 

The free energies AG^^, AG^^^ , AG^g^, and AGlj, were estimated using the thermody- 
namic integration (TI) method |15 | 119 1 ^18]. To minimize the numerical integration errors we 
employed the polynomial regression techniques to calculate free energy difference, instead of 
trapezoidal quadrature ^31J- Separate simulations were performed for changes in the Lennard- 
Jones with 21 values of the scaling parameter, A = 0.0, 0.05, 0.1 ...0.9, 0.95, and 1.0, and 
the electrostatics with 11 A values, A = 0.0, 0.024, 0.095, 0.206, 0.345, 0.5, 0.655, 0.794, 0.905, 
0.976, and 1.0. For simulations with only Lennard- Jones, all partial charges were set to zero 
and the soft-core scaling parameter was set to 0.5. Once the neutral atoms were fully grown 
in the solvent, the second simulations then computed the free energy associated with the elec- 
trostatics with a soft-core scaling parameter of 0.0. This was accomplished by increasing the 
partial charges from zero to their final values given by the forcefield. 

The free energy associated with the restraints, AG^^ was calculated using the Bennett 
acceptance ratio approach [32]. We performed 1.0 ns equilibrium simulation for the estradiol- 
receptor complex using each of the harmonic restraining potentials with force constants of 0, 25, 
40, 60, 90, 150, 200, 300, 450, 700, and 1000 kJ/mol/nm^ for distance, kJ/mol/rad for angle, 
and kJ/mol/rad^ for dihedral restraints. The first 0.5 ns of each simulation was discarded for 
equilibration and the remaining 0.5 ns was used to compute the free energy differences. No 
attempt was made to optimize the efficiency of the calculation since our primary objective was 
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Salmonid Whole 




Figure 4: Schematic showing the inferred evolutionary history of the Oncorhynchus mykiss 
estrogen receptors. Vertical bars mark inferred whole genome duplication events; short branches 
mark inferred duplications that were lost over time. 

to obtain accurate estimates of the restraining free energies. 
Evolutionary Analyses 

The following sequences were extracted from GenBank: AB037185, AF349412, A133920050, 
AY727528, AY775183, BD105560, AB190289, AJ487687, AY055725, AF061275, AF253505, 
AY520443, AJ242741, DQ009007, DQ248228, DQ177438, X89959, TNU7560, AY422089, AF298183, 
AF136979, AY074780, AB007453, AJ006039, AF253062, AY223902, ORZMER, AY917147, 
AF326201, AY305026, NM_180966, NM_174862, AB003356, AB070630, AB070901, AB083064, 
AB117930, AB190290, AF061269, AF136980, AF177465, AF185568, AF298181, AF298182, 
AF349413, AF349414, AF516874, AJ275911, AJ289883, AJ314602, AJ314603, AJ414566, AJ414567, 
AJ489523, AJ580050, AY074779, AY211021, AY211022, AY305027, AY307098, AY508959, 
AY566178, AY770578, AY917148, BC044349, BC086848, DQ177439, DQ248229, TNU75605. 
The first 30 are ERa sequences and the other 39 are EKP sequences, and the following analysis 
was done separately for these two subtypes. The codons were aligned based upon their aligned 
amino acid sequences, and these alignments were used to infer tree topologies using the neighbor 
joining method. Then the ligand binding domains were extracted from the alignments. PAML 
was used to test several codon-based likelihood models that allow for variable dN/dS ratios 
among lineages based upon the inferred phylogenies and the aligned ligand binding domains 
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Results and Discussion 



Experimental binding affinity results are not readily available for the four trout ERs due to 
the difficulty in isolating the different isoforms. Thus, to verify our methodology for estimating 
AGbind for rainbow trout ER-E2 we first performed simulations using human ER (PDB: IQKU) 
at 300 K and compared the binding affinities to the experimental results. Our computational 
estimates at 300 K are —63.2 kJ/mol for insertion (when interactions between E2 and its 
environment are turned on) and —59.6 kJ/mol for deletion (when interactions between E2 and 
its environment are turned off). Experimental binding affinity for human ER is —52.3 kJ/mol 
at 300 K [34J. Thus, our human ER binding affinity estimates are within about 10 kJ/mol of 
experiment which is within the expected error due to the atomic models [22]. The trout ER 
simulations followed exactly the same procedure as human, beginning with docking the E2 into 
the ER. It is important to note that our ER simulations were performed at 277 K to closely 
mimic the water temperature of rainbow trout natural habitat. 

Table [T] shows our binding affinity results from both insertion and deletion. Simulation re- 
sults from both deletion and insertion of electrostatics and Lennard-Jones interactions provide 
a rudimentary assessment of the accuracy of our calculations (note that there may be inaccu- 
racies in the atomic models but that is beyond the scope of this study). The fact that both 
insertion and deletion give very similar results strongly suggests that our simulations are well 
converged and that accurate estimates of binding affinities have been obtained. 

Our results in Table [l] show that the E2 binds preferentially to the ERal isoform of the a 
subtype that has been found in all salmonids. The other isoform ERa;2, which appears to have 
arisen during the recent salmonid whole genome duplication, shares 75.4% sequence identity 
with the ERofl and thus a large number of substitutions have accumulated since the initial 
duplication event. To infer the evolutionary pressures that led to this amount of divergence 
in both protein sequence and function, we examined the lineage specific differences in dN/dS 
ratios among the ERa sequences. We used an alignment of the codons in the ligand binding 
domain for all ERa sequences and a phylogeny inferred from the nucleotide sequence by the 
neighbor joining method (which did not differ significantly from the tree in [7]). PAML was 
used to calculate the log likelihood values and dN/dS ratios for each of five hypotheses: a single 
ratio for all branches, one ratio for all branches except the branch to the rainbow trout ERa2, 
separate ratios for the two ERa's from rainbow trout and the rest of the tree, separate ratios for 
the rainbow trout ERa2, all ERal from salmonids and the rest of the tree and the full model 
where every branch has its own ratio (see Table [2]). Using the Aikaike Information Criterion, 
the model with two ratios, one for the branch to the rainbow trout ERa2 and one for all other 
branches is the best fitting model. For this model, the dN/dS ratio for all other branches was 
0.09 whereas the ratio for the ERa2 branch was 0.30. In all tests, the dN/dS ratio for the 
ERa2 branch was about three times greater than the other salmonid branches. Therefore, the 
ERa2 ligand binding domain appears to be evolving under relaxed selection relative to the 
other salmonid ERal ligand binding domains, which is consistent with the decreased affinity of 
this domain for E2. It is also possible that ERa2 was evolving in a neutral fashion for a short 
time, but then developed a new function and is now undergoing stronger purifying selection. 
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Table 1: Estrogen receptor binding affinities for different isoforms obtained at 277 K. All results 
are in kJ/mol. The binding affinities AGbind were calculated using Eqn[T] Both insertion and 
deletion directions give very similar results which demonstrates that our simulations are well 
converged. 



(A) Insertion 





ERal 


ERa2 


ER/31 


ER/32 




-21.5 


-27.6 


-31.4 


-24.2 




-146.7 


-126.5 


-132.4 


-144.1 




48.9 


49.9 


52.4 


49.9 


AGi, + AGfj 


61.8 










-6.8 








^^bind 


-64.3 


-49.2 


-56.4 


-63.4 


(B) Deletion 




ERal 


ERa2 


ER/31 


ER/32 




-21.5 


-27.6 


-31.4 


-24.2 


AG2f, + AG|f 


-149.2 


-125.2 


-131.1 


-144.7 




48.9 


49.9 


52.4 


49.9 


AGi, + AGfj 


62.0 










-6.8 








^^bind 


-66.6 


-47.7 


-54.9 


-63.8 



This possibility could be explored further if more ERa salmonid gene sequences were made 
available. 

Our results show that both ER/3 isoforms bind similarly to E2, i.e., the difference between 
them in binding affinity is small compared to the difference between the ERa isoforms (see 
Table [T|. The two isoforms share only 57.6% sequence identity, having arisen prior to the 
Teleost radiation, and the difference in their binding affinity might be expected to be greater, 
given this large degree of divergence. We performed a similar analysis of the dN/dS ratio for 
these genes by testing the following models: one dN/dS ratio for the whole tree, a dN/dS ratio 
for each of isoform ER/31 and ER/32, dN/dS ratios for each of the two rainbow trout isoforms 
and for each isoform for all other fish and the full model where every branch has a different 
dN/dS ratio (Table [2]). The best fitting model for this comparison was the single dN/dS ratio 
(0.07) for the entire tree, indicating that both ER/3 isoforms are under the same level of purifying 
selection. This is also consistent with our results showing that these two ligand binding domains 
have similar affinity for E2. 

These nuclear ERs have a significant and ubiquitous distribution in the rainbow trout [71 [5] • 
The levels of transcription differ among the four genes with one isoform having higher transcript 
levels in most tissues than the other isoform. For the ERa isoforms, ERal has the higher 
transcript levels, and for the ER/3 isoforms, ER/32 has the highest transcript levels [7]. While 
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Table 2: Results of fitting evolutionary models for differences in dN/dS ratios. LBD and DBD 
indicate ligand and DNA binding domains, respectively. Oal, Oq2, 0/31 and 0/32 are the O. 
mykiss ERal, ERa2, ER/31 and ER/32 genes, respectively. Sal indicates all of the salmonid 
ERqI genes. /31 and /32 indicate ER/31 and ER/32 from all fish, respectively, np is the number 
of parameters in the model, InL is the log likelihood calculated by PAML, and AIC is the 
Akaike Information Criterion value [35]. Models labeled with an asterisk are the best fitting 
models based upon the AIC values. 



ERa-LBD 


np 


InL 


AIC 


Hq: Everyone is equal 


59 


-5645.4 


11409 


Hi: Oa2 / others 


60 


-5641.7 


11403 * 


H2: Oa2 / Oal / others 


61 


-5641.6 


11405 


H3: Oa2 7^ Sal 7^ others 


61 


-5641.5 


11405 


HfuII- Everyone is different 


115 


-5610.3 


11451 


ERa-DBD 


np 


InL 


AIC 


Hq: Everyone is equal 


59 


-1121.0 


2360 * 


Hi: Oa2 / others 


60 


-1120.8 


2362 


H2: Oa2 / Oal / others 


61 


-1120.7 


2363 


H3: Oa2 7^ Sal 7^ others 


61 


-1120.3 


2363 


HfuII- Everyone is different 


115 


-1101.8 


2434 


ER/3-LBD 


np 


InL 


AIC 


Hq: Everyone is equal 


77 


-7064.6 


14283 * 


Hi: /32 7^ /31 


78 


-7063.9 


14284 


H2: 0/32 / 0^1 7^ /31 7^ /32 


80 


-7062.1 


14284 


iJpuii: Everyone is different 


151 


-7005.9 


14314 



the correlation between reduced transcription levels and binding affinity is clear in the ERa 
isoforms, there seems to be no such correlation for the ER/3 isoforms. These two isoforms share 
similar binding affinity, and yet, ER/31 has much lower expression levels than ER/32 in juvenile 
rainbow trout. It is possible that both ERq2 and ER/31 have higher expression levels at other 
life stages [7J. Given the age of ER/31 and the equivalent levels of both E2 binding affinity and 
purifying selection compared with ER/32, this ER clearly continues to have an important role 
as an estrogen receptor. 

It is not as clear what ERa2's role is as an estrogen receptor. It's reduced affinity for E2, 
low transcript levels and evidence for relaxed selection suggests that this estrogen receptor may 
be undergoing subfunctionalization or neofunctionalization. One indication that ERa2 may be 
undergoing neofunctionalization is that the DNA binding domain of ERa2 does not have the 
degree of sequence variation that the ligand binding domain has. If the ERa2 was undergoing 
relaxed selection along it's entire length, the DNA binding domain would also show indications 
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of greater amino acid divergence (Table [2]). It appears that ERa2 is not losing its ability to 
bind to the canonical estrogen receptor element even though it is losing affinity for E2. This 
suggests that this gene may be undergoing neofunctionalization by binding to some other ligand 
than E2. 

Conclusions 

Using molecular dynamics simulations we estimated the binding affinities between the hor- 
mone 17/3-estradiol (E2) and different estrogen receptor (ER) isoforms in the rainbow trout, 
Oncorhynchus mykiss. Our results show that E2 binds preferentially to ERal over ERa2. A 
recent genome wide duplication event led to two functional ERa isozymes in O. mykiss. Our 
evolutionary and functional analyses along with Nagler's evaluation of transcription levels [7] 
suggest that the ligand binding domain of ERa2 has been or is currently evolving under re- 
laxed selection relative to ERal. Low sequence divergence of its highly conserved DNA binding 
domain suggests that ERa2 is likely undergoing neofunctionalization, in which it continues to 
recognize the same estrogen receptor element in the DNA but may be binding to a different 
ligand. For the ER/3 subtype both isoforms bind similarly to E2, in keeping with our evolution- 
ary analyses that both isoforms of this subtype are evolving under the same degree of purifying 
selection. 
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